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ABSTRACT 

The  small  perturbation,  two-dimensional  transonic  equation 
is  manipulated  with  a  separation- of -variables  approach  to 
obtain  two  ordinary,  nonlinear,  differential  equations. 
Numerical  integration  of  these  differential  equations  results 
in  new  transonic  boundary  surfaces  for  planar  external  flows. 
A  key  ingredient  in  these  solutions  is  the  identification  of 
dependence  of  two  integrations  constants,  a  and  0,  on  the 
parameter  (1-M^2)  .  The  anticipated  behavior  for  both  the  Mach 
number  and  the  pressure  coefficient  is  used  as  a  guide  in  the 
actual  selection  of  the  adjustable  constants  in  the  problem. 
The  physical  reality  of  our  boundary  surfaces  is  examined  by 
displaying  the  boundary  conditions  they  satisfy.  The  strictly 
sonic  flow  (Moo  =  1.0)  has  an  analytic  representation 
corresponding  to  a  divergent  surface  which  goes  supersonic. 
This  sonic  solution  is  compared  with  an  Euler-CFD  approach 
confirming  the  validity  of  our  results  over  the  region  where 
small  perturbations  apply.  Solutions  are  also  shown  for 
Moo  =  0  .  8  ,  0  .  9  ,  1 . 1 ,  and  1.2  .  These  results  are  consistent  with 
known  behavior  for  both  subsonic  and  supersonic  external  flow. 
Since  the  results  of  this  work  yield  actual  transonic 
contours,  we  can  examine  shockless  surfaces  for  design 
applications.  The  possibility  of  starting  with  transonic 
surface  is  of  interest  to  present  day  CFD  approach.  Finally  an 
entire  transonic  upper  surface  is  presented  for  Moc  =  0.8,  by 
patching  a  subsonic  Mach  number,  which  reaches  a  plateau  at 
M=1.0,  with  a  sonic  flow.  This  patching  requires  the  careful 
interpretation  of  a  nondimensional  reference  length,  called 
Y0,  which  is  a  function  of  M„. 
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I .    INTRODUCTION 

The  governing  nonlinear  equations,  together  with  realistic 
boundary  conditions  have  always  been  the  main  difficulty  in 
transonic  aerodynamic  analysis,  [Ref .  1] .  While  the  hodograph 
transformations  approach  linearizes  these  equations,  the 
transformations  are  confined  to  a  limited  number  of  transonic 
flows  due  to  complicated  boundary  conditions.  The  small 
perturbation  approach  simplifies  the  full  potential 
Equations;  furthermore,  like  other  predictive  methods,  it 
succeeds  in  retaining  the  great  sensitivity  of  transonic  flows 
within  small  perturbations.  This  paper  expands  on  the  exact 
solution  of  the  small  perturbation,  nonlinear,  two 
dimensional,  transonic  equation,  using  the  guidelines  and 
tools  given  by  Biblarz  0.,  [Refs.  2  and  3]  .  Starting  with  a 
separation-of -variables  approach,  Chap.  Ill  reveals  two 
nonlinear,  ordinary,  differential  equations  which  lead  into 
two  exact  solutions.  Chaps.  IV  and  V  describe  the  numerical 
integration  of  these  implicit  solutions,  which  finally  yield 
the  boundary  surfaces.  Later,  these  boundary  surfaces  are 
found  to  satisfy  the  boundary  conditions  for  a  two  dimensional 
surfaces.  These  findings  are  then  compared  with  a  finite 
difference  solution  of  the  Euler  Equations  for  the  sonic 
solution.  The  usefulness  of  the  shockless,  transonic,  boundary 
surfaces   appears   once   transformed   to   dimensional   body 


surfaces;  "patching"   (i.e.,  translating),  the  dimensional, 
subsonic,  body  surface  of  Cp=0  at  the  upstream  inflow,  with 
the  dimensional,   sonic  body  surface,   enables  us  to  form 
complete  transonic  upper  surfaces  of  interest  for  design. 


II.    THE  TRANSONIC  EQUATION 

The   small   perturbation   nonlinear,   two   dimensional 
Transonic  Equation 

(1-M„2)$   +  $   =  M~{t+1)  $   $  (1) 

is  described  by  the  velocity  potentials 

u  =  <&x  v  =   $y  (2) 

Multiplying  eqn.   1  by  1/Ua  enables  us  to  transform  the 
Transonic  Equation  to 

(l-M2)^  +  ^  =  ^(y  +  DiM)^  <3> 

where  the  nondimensional  velocity  potentials  are 


<t>  =  —  <t>  =  —  (4) 


The  Meyer  solution  for  the  de- Laval  nozzle  defines  a 
velocity  potential 


<t>  (x,y)    =  b 


—  +  b(y+l)M2-&£  +  b2  {y+\)2Mt^- 


(5) 


which  satisfies  the  nondimensional  form  of  equation  (3)  for 
Moo  =  1.0  .  The  Meyer  solution  belongs  to  a  general  class  of 
sonic  flows  as  reported  by  Guderley  and  Yoshihara,  [Ref.  4]. 


The  transonic  equation  (1)  may  be  modified  to  a  more 

useful  form  if  multiplied  by  M„2  (7  +  1)  /U&, ,  instead  of  just  1/U- 

„,    to  give 

(l-MJ)^  +  4>yy  =  (J>x4>xx  (6) 

Here  the  modified  velocity  potentials  are 

4>x  =  Mhy  +  D-jy-  4>y  =  Afi(Y+D-^     (7) 

For  a  given  problem,  Mw  and  7  are  constants,  and  eqn.  6  is 
really  less  cluttered.  However,  notice  that  Ma2  (7  +  1) /Ua  is 
already  given  in  the  derivatives  of  <p  in  eqns .  7,  and 
therefore  appears  in  the  expressions  of  Cp  and  the  boundary 
surfaces  as  will  be  shown  later.  This  form  of  the  modified 
transonic  equation  (6)  will  be  examined  in  the  next  chapters. 


III.     THE  GOVERNING  DIFFERENTIAL  EQUATIONS 

An  exact  solution  to  the  modified  transonic  equation  6  has 
been  given  by  Biblarz,  [Ref.  2],  starting  from 

4>(x,y)    =<bs(x,y)    +  (l-Afi)x  (8) 

where  <£s(x,y)  refers  to  the  sonic,  small  perturbation  solution 
(M00  =  1.0)  .  In  other  words,  "any  flow  which  satisfies  the  sonic 
equation  also  satisfies  the  full  transonic  range  with  similar 
boundary  conditions,"  [Ref.  2]. 

An  exact  solution  may  be  obtained  by  using  the  separation- 
of -variables  approach  with  the  potential  function  $(x,y)  of 
the  form 

(J)(x,y)  =£(x)Ti(y)  +  (l-Mi)x  (9) 

Substituting  the  above  <p  function  in  the  modified 
transonic  equation  6  results  in  two  ordinary,  second  order, 
nonlinear  differential  equations 

ip.^11  -  \l   =  o  (10) 

dx  dx2 


^     ~    ATI2  =  0  (11) 

dyl 


where  X  is  a  separation  constant. 


A  simple  solution  to  the  first  differential  eqn.   10  is 
obtained  by  multiplying  both  sides  by  d£/dx, 


dx\  dx  ax2 )        dx 


(12) 


or 


d 
dx 


l/^i\3  -  1^ 


3  dx 


=   0 


(13) 


Thus 


&  - 

dx        \ 


2X^2 


+   a 


(14) 


dx   = 


dt, 


\ 


3XV 


+   a 


(15) 


and 


x  -  xn 


■/■ 


tfS 


3A 


S2  +  a 


(16) 


where  a  and  x0  are  integration  constants. 

A  similar  procedure  can  be  performed  to  solve  the  second 
differential  eqn.  11  yielding 


tfyidy2J   dy 


(17) 


or 


d 
dy 


2\dy/   3  ' 


=  0 


(18) 


t  -^ 


3     ' 


(19) 


dy 


dn 


Nfi3  +  P 


(20) 


and 


y  -  y0  = 


=/• 


dn 


\ 


^n3  +  p 


(21) 


where  jS  and  y0  are  integration  constants. 

Eqns .  16  and  21  are  exact  solutions  to  the  differential 
eqns .  10  and  11.  In  fact,  boundary  conditions  are  needed  to 
evaluate  the  integration  constants  a,  /5,  Xq,  y0  as  well  as  X. 
An  important  step  in  this  procedure  is  the  recognition  of  ce 
and  /?  as  functions  of  (1-Mj)  (  [Ref.  2],  i.e., 

a   =   ±q  (1-Mi)2  (22) 


p  =  ±c2  (1-Mi) 


(23) 


The  upper  positive  sign  and  the  lower  negative  sign  in  the 

above  equations  will  be  shown  to  belong  to  M^sl.O  and  Moo<1.0 

respectively,  namely, 

asO  &  0sO     for      M.&1. 0 
a<0  &  /5<0      for      M(X<1.0 

See  Appendix  B. 


The  constants,  a,    /?  and  X  can  be  shown  to  be  related  by  the 
expression 


Sf    -   |<1-«.V 


Since  C,  ,  C2  and  X  are  positive  constants 


1  2    3 


(24) 


(25) 


The  usefulness  of  these  definitions  will  be  outlined  in 
more  detail  once  we  use  the  exact  solutions  from  eqns .  16  and 
21,  while  verifying  the  correct  boundary  conditions  for  the 
pressure  coefficient  Cp. 

Upon  inserting  eqns.  22  and  23  in  eqns.  16  and  21  we  have 


x  -  xr 


■/■ 


di 


\ 


i±l2  ±  q(i-M»2)2 


(26) 


y  -  y 


.-*/■ 


dr| 


2\ 
NT 


ti3  ±  c2{i-mI) 


(27) 


Factoring  out  Cj  and  C2  we  obtain 


X     -     Xr, 


-  O  / 


c?S 

\ 

3*  £2   ±    (1-Mi)2 

(28) 


y  -  y0 


±  c2   2 


/■ 


dTl 


X  X]3   ±  (1-M2) 


N  3C: 


(29) 


Introducing  new  dependent  variables 


5  -  5 


3X 


U 


2C 


i/ 


(30) 


ii  =  n 


(31) 


the  above  simply  becomes 


x  -  x, 


\  /3A\-i    r  ^ 


•  ■  ^  (f  r  fe 


^2  ±  (1-mH) 


(32) 


y  -  y0  = 


■ii2X- 


cffj 


±  (i-m: 


(33) 


Finally,  if  we  define  new  independent  variables 


X*    (x  -*b>Ci~*(-^ 


(34) 


r  =  (y  -y0)c^(-^" 


(35) 


eqns .  32  and  33  then  follow  as 

dl 


*-/■ 


j¥   ±    U-Ml) 


(36) 


y-*//3   J    2  <37) 

This  set  will  be  numerically  integrated  to  obtain 
transonic  surfaces  with  their  corresponding  Cp  and  Mach  number 
profiles . 
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IV.    NUMERICAL  INTEGRATION 

We  can  rewrite  eqns .  36  and  37  in  a  more  general  form  and 
specify  the  two  integral  limits 

±1 

X  =   /  1 ~ (38) 

o  Jz2  ±  (i-Mi)2 


Y       ±         f 


dW__ 

2 


?d-w.2)  3  v 


As  indicated  before,  the  upper  and  lower  sign  within  the 
integrand  root  and  the  integral  limits  is  used  for  M^sl.O  and 
Moo<1.0  respectively.  In  addition,  our  flow  of  interest  is  only 
the  first  quadrant  in  X  and  Y,  which  represents  physical 
external  flows.  To  show  this  clearly,  for  the  subsonic  case 
M^<1 . 0 ,  eqns.  38  and  39  reduce  to  the  form 


X  =  f-3 ^ OzU-(l-MH)  (40) 

0   y/z2  "  (1-M„2)2 


y=-    f      ^ f\*(i-*d)^ 

JW3    -    il-Ml) 


(41) 


(l-AC) 


2\     3 


11 


and    for  Mw2l.O    we   have 


*-/l — l*°  (42) 

o      sjz2    +    {1-Ml)2 


(1-M.)    3 


dW 

^3     +      (1_MJ) 


r=        f        ^ fj.-d-Mi) 


(43) 


However,   we  will  frequently  retain  the  more  general  "±" 
notation  in  our  discussions  further  on. 

Before  we  head  to  numerical  solutions,  we  should  indicate 
that  sonic  flow  (Moe  =  1.0)  can  be  given  in  explicit  form  through 
eqns .  40,41,42,43,  i.e.,  the  resulting  equations  for  sonic 


flow  are 


I 


*• 


dZ 


(44) 


'o  W* 


dW 


Y=-f^L  (45) 


Solving  equations  44  and  45  we  obtain  our  sonic  solution 

I   =(\A2  <46) 


-f 
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Having  defined  |,^,X  and  Y  in  eqns .  30,31  and  eqns .  34,35, 
we  can  express  the  sonic  solution  rewriting  eqns.  46  and  47  as 

5  =  ^(^'(yl  (48) 


(y-vXl^  <49) 


Providing  the  sonic  perturbation  potential  from  eqn.  9 

4>5  =  5n  (so) 

and  inserting  eqns.  48  and  49  we  have,  [Ref.  2] 

1  (X-X0)3 

3  (y-y0)2 

Next,  eqns.  40  to  43  are  numerically  integrated  and 
plotted  in  Figs.  1  and  2  for  various  MB,  (e.g., 
0.8,0.9,1.0,1.1,1.2)  .  Plotting  rj  vs.  Y  for  Y>0  in  Fig.  2  is  an 
appropriate  choice  for  investigating  physical  external  flows 
as  stated  earlier. 

Our  numerical  solution  has  been  based  upon  the  Newton- 
Cotes  method  of  order  4,  [Ref.  5] ,  which  may  be  applied  to  an 
even  number  of  subintervals ,  each  described  as, 

J  fix)  dx=V^Li\7  f(a)  +32  ^p)+12  ^(^)+32  f(Al^-)+7  f(b) 

(52) 

(See  appendix  A) . 
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Numerical  solutions  which  express  f  and  rj ,  as  has  been 

noted,  can  be  combined  with  eqns .  30,31  and  eqns .  34,35  and 

substituted  in  eqn.  9  to  give  a  complete  solution  for  the 

perturbation  potential  in  terms  of  the  positive  constants  X  or 

(C1;  C2)  ,  and  (Xq,  y0)  . 
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V.    BOUNDARY  SURFACES 

For  an  inviscid  flow,  the  condition  to  be  applied  at  the 
surface  of  a  solid  boundary  is  that  the  direction  of  the  flow 
velocity  vector  is  tangent  to  the  solid  surface,  [Ref .  6] .  In 
terms  of  perturbation  velocities  this  boundary  condition 
becomes 

\dx '/  surface  U~ 

The  modified  velocity  perturbation  potential  as  pointed  out 

earlier  is 

<j)y  =  M„2(y+1)-^  (54) 


thus 

v  $, 


u»        Mi(Y+D 
We  can  show  by  substitution  that 

dy\  _  <t> 


(55) 


y—  (56) 

CtX  I  surface  M^ij+l) 


Now,  deriving  the  perturbed  potential  function  in  eqn.  9 
with  respect  to  y  gives 


dr\ 

dy 


4>y  =  e-?L.  <57) 
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Eqn.    19    for   y>0    becomes 


lAn3   +   p  (58) 


~dy         \~2 

This  relation  can  be  rearranged  using  r\   and  /5  in  eqns .  31  and 
23  to  yield 

41  -  -C?J  ^   ±  d-Mi)  <59> 

dy       v 

Upon  using  this  equation  and  the  relation  of  eqn.  30  to  the 
derived  velocity  potential  in  eqn.  57  results  in 

4>y  =  -I  (-!£]"*  Q^/fp  ±  d-Mi)  (60) 

The  last  equation  can  be  rewritten  with  the  benefit  of  the 
related  constants  C, ,  C2  and  X  in  eqn.  25,  hence 


(j)y  =  -2ty  fj3  ±  il-Ml)  (61) 


Once  we  have  found  the  form  of  <f>  ,    the  expression  for  the 
surface  boundary  condition  is  given  as 

(#)      --T-T^ W^3  ±  <1"^>         (62) 

\  ^/surface  3    M^Y+D 

However,  an  important  result  to  observe  is  that  £  and  rj  in  the 
R.H.S.  of  this  equation  are  given  implicitly  as  a  function  of 
X  and  Y  respectively  in  our  earlier  numerical  integration 
solution,  while  the  L.H.S.  is  expressed  by  x  and  y.  To  get 


16 


compatibility,  we  then  transform  the  L.H.S.  of  eqn.  62  using 

the  chain  rule  ,     ,    ,  .  , 

dy   .  dY    dX/dx  (63) 

dx   "  dX    dY/dy 

Deriving  eqn.  34  with  respect  to  x  and  eqn.  35  with  respect  to 

y  gives  x  i 

d£  =   r  ""6/lA\'2  (64) 

dx  x  2 


dy       C2  I  3  ' 


therefore 

^  =  24?  (66) 

dx        2  dX 

Equation  62  can  be  rewritten  using  the  above  transformation  as 

<**  '3'  Mi(y+1) 


or 

dY 


y/   fj3  ±  (1-MM2)     l3/  ^(y+D 


E  <*  (68) 


We  are  now  in  a  position  to  compute  the  exact  boundary 

surfaces  out  of  eqn.  68.  Starting  with  the  L.H.S.  we  then 

define 

E(Y)  (69) 


yj   fj3    ±     {1-Ml ) 
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Having  defined  eqn.  69  we  can  plot  E  (Y)  vs.  Y  in  Fig.  3  for 
various  Mw,  (e.g.,  0.8,0.9,1.0,1.1,1.2),  keeping  in  mind  our 

"±"  notation  for  the  supersonic  (upper  sign)  and  subsonic 

(lower  sign) . 

The  asymptotic  value  of  E (Y)  occurs  at  Y0  where 

fj^(l-Mi)^  (70) 

In  other  words,  the  asymptotic  value  can  be  expressed  as 

Y0=Y0(Moo),   or   as   function   of   (1-Ma2)  ,   which   has   been 

numerically  found  to  be 

2.4 

Y     s  

0  i  (71) 

|1-M„2|6 

Y0  vs.  M„  is  sketched  in  Fig.  4  to  demonstrate  an  interesting 
symmetry  with  respect  to  Moo  =  1.0  .  Here  Y0=oo  but  note  that  for 
practical  reasons  we  can  adopt  finite  values  of  Y0  (perhaps  as 
low  as  7  to  20)  to  represent  conditions  sufficiently  close  to 
Moa=1.0  . 

Defining  the  integral  of  the  R.H.S.  of  eqn.  68  as 


K(X)     =   "(If— -i f  I   dX  (72) 

\3/  AL2(V+1)  Jo 


and  redoing  numerical  integration  allow  us  to  plot  K(X)  vs.  X 
for  7=1.4  in  Fig.  5.  It  is  important  to  realize  that  K(X)>0 
for  M^cl.O  and  K(X)<0  for  M^sl.O  ,  before  solving  eqn.  68. 
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Therefore,  integrating  both  sides  of  eqn.  68  yields  for  K(X)>0 

I'd 

(E(Y)dY  =  K(X)  (73) 


and  for  K(X)<0 


y 


[E(Y)dY   =  K(X)  (74) 

Once  we  specify  X  for  K(X)  we  can  immediately  solve  the  L.H.S. 
integrals  of  eqns .  73  and  74  for  Y,  expressing 

(Y)  surface    =  *<*)  <75) 

namely,  the  required  boundary  surface  for  any  chosen  M„. 

We  can  rewrite  eqns.  73  and  74  in  the  more  complete  form 
for  IVL<1 .  0  as 


{  yj   fj3  -  (1-Mi)     ^3/  Mliy+1)    { 


dX  (76) 


and  for  M^sl . 0 

Y 

dY 


(77) 


Eqns.  76  and  77  are  equivalent  to  eqns.  73  and  74  which 
finally  yield  the  boundary  surfaces.  These  are  shown 
nondimensionaly  in  Figs.  6  and  7  for  Moo  =  0  .  8  ,  0  .  9  ,  0  .  97  and  in 
Fig.  8  for  M0o  =  l .  0  ,  1 . 1,  1 .  2  . 
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The  sonic  surface  is  most  easily  verified  by  direct 
analytical  expansion  of  eqn.  77,  using  the  sonic  relations 
stated  in  eqns .  46  and  47. 

/(ip*  ■-  [if  w/(l^ 

Solving  the  above  integrals  on  both  sides  gives 

r4  -  r4  =  -   i   (l\5x4  (79) 

Recall  that  Y0  -*  oo  for  the  sonic  case.  However,  choosing  for 
example  Y0=20  gives  an  asymptotic  value  of  17  =  0.01,  (eqn.  47), 
and  E(Y0)=1000,  (eqn.  69),  for  our  explicit  sonic  solution. 
(See  also  Figs.  2  and  3) .  Therefore,  dividing  by  a  finite  Y0 
enables  us  to  obtain  the  nondimensional  sonic  boundary  surface 
as  a  function  of  y    . 


/2\5(  A")4  (80) 


r0   N1"  (Y+1)  ^3/  \Yo 


This  will  be  shown  to  be  a  divergent  surface  which  increases 
the  Mach  number  to  supersonic  conditions. 
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VI.    PRESSURE  COEFFICIENT 


For   a   flow   with   small   perturbations   the   pressure 
coefficient  is  given  by,  [Ref.  8] 


Cp  =  ~ 


^  +  {1-M^)(JL)\(^V' 


(81) 


A  first  order  approximation  for  the  linearized  pressure 
coefficient  in  a  two  dimensional  planar  flow  is 

2u 


C*--u 


(82) 


Recall  the  modified  velocity  potential 


<bx  =  al2(y+d 


u 


(83) 


Thus 


u 


4>, 


u~        Ad(Y+D 


(84) 


and  from  eqn.  82  we  see  that 


Cp=  ~ 


2(b 


m:(y+i 


(85) 


Now,   deriving  the  potential  function  in  eqn.   9  with 
respect  to  x  we  have 

dx 
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Rewriting  eqn.  14  using  the  relation  for  a   in  eqn.  22  as 


dx   \ 


3\ 


V  ±   q(i-Mi) 


2\  2 


(87) 


By  taking  advantage  of  £  as  stated  in  eqn.  31  and  factoring 
out  C,,  eqn.  87  becomes 


-^  =  cj   V  V  ±  d-^-2) 


(88) 


Substituting  eqn.  88  and  the  relation  for  rj    in  eqn.  31  to  the 
velocity  potential  </>x  in  eqn.  86  yields 


4>*  =  fl 


(i^)*  7FT1 


i-Mi)"  +  (i-m!: 


(89) 


Once  again,  using  the  benefit  of  the  related  constants  C, ,  C2 
and  X  in  eqn.  25,  eqn.  89  becomes 


4>x  =  fj  V  V  ±    il-Ml)2   +  (i-nd) 


therefore  eqn.  85  for  C  follows  as 


(90) 


C?=     —2 


m:(y+d 


fj  J  l2   ±    (1-M„2)2  +  (l-Mi) 


(91) 


We  must  anticipate  the  above  expression  for  Cp  to  satisfy  the 
boundary  conditions  on  external  boundary  surfaces  as  will  be 
explored  further  on. 
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VII.     THE  LOCAL  MACH  NUMBER 

The  velocity  vector  V  is  given  by  the  undisturbed  uniform 
velocity  U0  and  the  perturbation  velocities  for  the  2D  flow  as 

V  =  I(Um  +  u)    +  jv  (92) 

Thus,  we  may  obtain  the  local  Mach  number  M  in  terms  of  the 
above  velocities  and  the  local  speed  of  sound  "a". 

\V\2          {UM  +  u)2  +  v2  .__. 

M2  =     I     I      = (93) 

a2  a2 

Neglecting  higher  powers  of  the  perturbation  velocities  gives 

/ 


Ul 


M 


2  _ 


i  +  M] 

U  (94) 


a2 


In  addition,  the  ratio  of  the  local  speed  of  sound  to  the 
undisturbed  uniform  speed  of  sound  becomes 


a2        _T_ 
ai  '  T. 


or 

*1   =  Tq/T" 
a2     '  TQ/T 


were  T0  designates  the  stagnation  temperature. 


(95) 


(96) 
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Eqn.  96  can  be  replaced  using  the  well  known  relation 


2 


1  +  1-Am2 
2 


(97) 


Substituting  the  local  speed  of  sound  from  eqn.  97  to  eqn.  94 

results  as 

1  +  -=± 


AT  =  AC 


'     Y-l   2N 
2 

1  +  ^Af2 
2 


(98) 


Having  stated  C  in  eqn.  82,  eqn.  98  is  reduced  to  the  form 


M< 


m:(i-cp) 


i  +  iz1m*      i  +  Jtzi 


m. 


(99) 


Rearranging  eqn.  99  and  factoring  out  M2  yields  the  final 
expression  for  the  local  Mach  number 


AC(1-C_) 

Mz   = 


1+  Y_iM„2Cn 
2     p 


(100) 


Note  that  M=M00  when  Cp=0  as  it  should. 

This  important  relation,  combined  with  the  corresponding 
eqn.  91  for  Cp;  will  be  used  when  discussing  the  boundary- 
conditions  in  the  next  chapter. 
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VIII.    BOUNDARY  CONDITIONS 

Rather  than  providing  the  boundary  conditions  for  the 
governing  differential  eqns .  10  and  11  to  evaluate  the 
constants,  we  imposed  a  solution  by  defining  a  and  (3  in  eqns. 
22  and  23  respectively.  That  is  to  say,  boundary  conditions 
have  been  already  implied  in  the  definitions  of  ce  and  13 .  This 
procedure  led  us  to  reveal  the  boundary  surfaces  upon  using 
the  tangency  condition  and  compute  the  pressure  coefficient  Cp 
on  the  body  surface,  (eqn.  91) . 

Thus,  we  arrive  at  the  necessity  to  verify  the  validity  of 
the  expression  for  Cp  for  the  boundary  surface,  and  verify  the 
following  requirements: 

•  At  the  upstream  end  (X/Y0=0)  ,  we  must  require  that  Cp=0 
ensuring  the  unperturbed  Mach  number  of  the  transonic 
inflow  to  equal  M0. 

•  For  the  subsonic  inflow  at  the  downstream  end  (X  @  Y/Y0=l) 
we  must  require  that  the  local  Mach  number  saturates 
towards  M=1.0  .  This  is  because  of  the  well  known  relation 
in  gas  dynamics  for  a  convergent  flow. 

We  first  restrict  our  attention  to  the  first  requirement, 

for  Mwal.O  at  X/Y0=0 

i  =  r  =  o  doi) 

At   Y/Y0=l    eqn.    70   becomes 

f|  -iT  =  -U-Mi)^  (102) 
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Inserting  the  above  terms  of  £  and  rj    in  eqns .  90  and  91  gives 

<|>x  =  -(i-m„2)3  ^{i-mI)2  +  (1-M„2)   =  0  (103) 

and 

cp  @  (x/y0  =  o  ,  r/y0  =  i)  =  o  (104) 

Therefore,    from  eqn.    100 

M®    (X/Y0   =  0    ,    Y/Y0   =  1)    =  Mm  (105) 

Hence,  eqns.  104  and  105  comply  with  the  first  requirement. 
The  complete  behavior  of  the  pressure  coefficient  and  the 
local  Mach  number  are  plotted  in  Figs.  9  and  10  for  Moo  =  1.0  and 
M00  =  l.l.  In  addition,  it  is  important  to  note  that  the  boundary- 
surfaces  ought  to  be  truncated  at  M^sl.2,  so  as  not  to  exceed 
the  small  perturbation  assumptions. 
Now,  for  Me, <  1.0  at  X/Y0=0 

I   =  I*  =   o  (106) 


However,    at   Y/Y0=0 

n  -*  fj 


(107) 


Thus,    eqn.    91    results   as 

Cp   ®(X/Y0   =  0    ,    Y/Y0   =  0)    -  *  (108) 

We  immediately  conclude  that  the  above  coordinate  (0,0)  can 
not  serve  as  the  upstream  start  for  the  subsonic  inflow. 
Instead,  we  must  localize  a  new  upstream  end  coordinate  along 
the  boundary  surface  whose  |*  and  rj*   values  satisfy  eqn.  91 
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for  the  equality  Cp=0  .  Since  every  coordinate  on  the  boundary 
surface  (X/Y0,Y/Y0)  is  a  function  of  |  and  rj ,  (eqn.  76)  ,  we  can 
numerically  scan  the  locus  of  all  the  coordinates,  i.e.,  the 
boundary  surface,  for  such  values  of  \*  and  rj*  that  will 
satisfy  eqn.  91.  For  example,  scanning  the  results  for  the 
subsonic  boundary  surface  of  1^=0.8  reveals;  £  =  £*=- 0  .  35763  and 
rj  =  rj,=  3.  01393  at  the  coordinate  (X/Y0=0  .  31366  ,Y/Y0=0  .40616)  . 
Inserting  the  above  values  of  £*  and  rj*  in  eqn.  91  yields 
Cp=0.0003  .  Hence,  by  relocating  a  new  upstream  coordinate  for 
the  subsonic  inflow  we  are  then  able  to  satisfy  our  first 
requirement.  (This  is  shown  later  in  Figs.  11  and  12) . 

We  now  consider  the  second  requirement  for  the  subsonic 
outflow  at  the  downstream  end  where 

I   -  T  =  -d-  Ml)  (109) 

fj  @(y/y0  =  i)  -  fj*  =  (i-mZ)1  {110) 

or  equivalently  dY/dX  -»  0 
Substituting  these  terms  for  |  and  rj    in  eqns .  90  and  91  gives 

(j>x  =  (i-mZ)  (in) 


and 

-2(1 -Mi) 
Cp  =  — (112) 

Mm2(y+D 
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Thus,  we  can  express  the  local  Mach  number  using  eqn.  100 


Ml 


'   2(1-Mj)  y 
1  + 


2     A*f(Y+l) 


Upon  expanding  eqn.  113  and  factoring  out  Ma2  we  obtain 


I2 
M2   =  '"  "  "  "  (114) 


AC(Y-D+2 

Afi(Y-l)  +  2 


or  equivalently 


M  =  1  (115) 


Hence,  we  have  verified  the  second  requirement:  The  local 
Mach  number  saturates  at  the  downstream  end  to  the  value  M=1.0 
for  any  subsonic  inflow  Moo<1.0  .  The  pressure  coefficient  and 
the  local  Mach  number  behavior  along  the  boundary  surface  are 
plotted  in  Figs.  11  and  12  for  Mw  =  0.8  and  Moo  =  0.9  respectively, 
showing  clearly  both  the  upstream  and  downstream  boundary 
conditions . 

Knowing  that  boundary  conditions  require  the  gradient  of 
the  potential  function  0  to  vanish  far  away  from  the  body 
gives  rise  to  an  important  question:  How  is  the  flow  pattern 
along  the  boundary  surface  propagated  into  the  flow  field?  We 
treat  this  question  by  considering  a  computationally  efficient 
flow  solver  developed  by  J.  A.  Ekaterinaris,  [Ref.  7].  This 
CFD  solver  consists  of  the  inviscid  Euler  equations  and  uses 
fully  implicit  two  dimensional  Crank-Nicholson  method  with  a 
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central  difference  AD  I  scheme,  and  fourth  and  second  order 
numerical  dissipation,  [Ref .  9] .  This  solver  has  been  modified 
here  to  solve  the  inviscid  flow  around  our  earlier  boundary 
surface  solution  for  the  chosen  sonic  case  Moo  =  1.0  .  By 
creating  a  grid  mesh  which  consists  of  a  boundary  surface 
edge,  we  are  able  to  map  the  flow  pattern.  Upon  observing  the 
constant  Mach  lines  as  graphed  in  Fig.  14,  we  deduce  that  the 
flow  is  shockless  through  the  whole  transonic  range;  a  shock 
is  formed  further  downstream  at  the  local  Mach  number  M=1.9  . 
These  results  are  in  complete  accordance  with  the  inviscid 
small  perturbation  transonic  solution,  (see  Fig.  9).  In 
addition,  each  vector  on  Fig.  13  describes  the  resultant  flow 
direction  at  a  specified  grid  point,  (not  the  streamlines  ! ) , 
colored  by  Mach  levels.  Thereby,  we  may  evidence  the  far- field 
boundary  conditions 

V  =  U„  u  =  v  =   0  (116) 

In  other  words,  the  flow  pattern  is  propagated  out  from  the 
boundary  surface  gradually  towards  satisfying  the  far- field 
boundary  conditions  far  away  from  the  body. 

It  is  of  special  interest  to  inspect  the  dimensional 
boundary  surface  for  possible  design  applications;  For 
example,  the  dimensional,  subsonic,  boundary  surface  of 
Moo  =  0.8  is  obtained  by  multiplying  the  nondimensional 
coordinates  (X/Y0r08),  Y/Y0(0g))  by  Y0(0.8)=2.84  ,  i.e.,  the 
dimensional,  downstream  end  coordinate  results  as  (X,,2.84), 
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as  plotted  in  Fig.  15a.  Similarly,  multiplying  the 
nondimensional ,  sonic  coordinates  (X/Y0(1 0),  Y/Y0(1 0))  by  Y0(10)=14, 
yields  the  dimensional,  sonic,  boundary  surface,  where  the 
dimensional,  upstream,  starting  coordinate  becomes  (0,14). 
Now,  translating  the  dimensional,  sonic,  boundary  surface, 
(Fig.  15b),  by  (AX, AY)  =  (X, ,  -  (14-2 . 84) )  ,  to  match  the 
downstream  end  coordinate  of  the  above  dimensional,  subsonic, 
boundary  surface,  "patches"  both  dimensional  boundary 
surfaces,  to  form  a  2D  planar  upper  surface  as  shown  in  Fig. 
15.  We  may  use  this  technique  to  "patch"  any  subsonic  boundary 
surface  in  order  to  form  our  transonic  upper  surface  of 
interest . 
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IX.    CONCLUSIONS 

We  have  shown  that  an  exact  solution  for  the  transonic 
equation  satisfies  the  complete  behavior  of  the  pressure 
coefficient,  and  the  local  Mach  along  a  2D-surface  for  planar, 
external  flow.  It  is  of  special  interest  to  inspect  the 
dimensional  boundary  surface  for  possible  design  applications; 
For  example,  the  dimensional,  subsonic,  and  sonic  boundary 
surfaces  are  obtained  by  multiplying  the  nondimensional 
coordinates  (X/Y0,Y/Y0),  by  the  corresponding  Y0.  Translation 
of  the  dimensional,  sonic,  boundary  surface,  to  match  the 
downstream  end  coordinate  of  the  above  dimensional,  subsonic, 
boundary  surface,  "patches"  both  dimensional  boundary 
surfaces,  to  form  a  2D  planar  upper  surface  as  shown  in  Fig. 
15.  This  technique  may  serve  as  a  powerful  tool  for  the 
implementation  of  supercritical,  shockless,  planar  surface 
designs,  because  of  the  saturation  of  rj  with  Y  (Fig.  2),  we 
can  also  describe  certain  internal  flows  such  as  planar 
converging,  diverging  tunnel  walls.  We  may  also  extend  these 
results  to  axisymetric  transonic  flow. 
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APPENDIX   A 
COMPUTER  PROGRAMS 

c*****  *************************  ***************************** 

*    THIS  FORTRAN  PROGRAM  USES  NUMERICAL  INTEGRATION  BASED  ON 
THE  NEWTON-COTES  FORMULA  TO  SOLVE  NUMERICALLY  EQNS .  40  AND  42. 

THE  OUTPUT  IS  IN  THE  FORM  OF  ^  vs.  X 
************************************************************ 

program  KSI 

PARAMETER  (N=5) 

IMPLICIT  REAL*8 (A-H,0-Z) 

IMPLICIT  REAL*8(M) 

IMPLICIT  REAL* 8 (K) 

INTEGER  INDEX, I 

DIMENSION  X(N) ,XDOT(N) ,SAVEX(N) , SAVED (N) ,KSI (N) 

DIMENSION  MACH (N) , XI (N) , X2 (N) 

MACH(l) =0.8D0 
MACH (2) =0.9D0 
MACH(3) =1.0DO 
MACH(4)=1.1D0 
MACH(5) =1.2D0 

KSI  (1)=0.0D0 
KSI(2)=0.0D0 
KSI (3) =1.0D-3 
KSI(4)=0.0D0 
KSI  (5) =0.0D0 

£*********************************************************** 

C  MAIN  PROGRAM  * 

c****  ********************************************  *********** 

DO  100  1-1,5 
IFQ.LT.3)  THEN 
H=-0.001D0 
GOTO  5 
END  IF 
H=0.01D0 
5    IF(I.LT.3)  THEN 

TF=- (l.ODO-MACH(I) **2) 
GOTO  10 
END  IF 
TF=120.0D0 
10   IF(I.LT.3.AND.KSI(I) .GT.TF.OR. 

I. GE. 3. AND. KSI (I) .LE.TF)  THEN 
30     do  50  index=l,5 

CALL  DERIVat ive ( KS I , N , I , XI , X2 , XDOT , MACH ) 
CALL  NEWTONCOTES (KSI, XDOT , X , S AVEX , SAVED , H , N , INDEX , I ) 
50    end  do 
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WRITE(10+i, *)  X(I),KSI(I) 
GOTO  10 
END  IF 
100  END  DO 
STOP 
END 
c*********************************************************** 

c  SUBROUTINE  DERIVative  * 

c**  *********  ************************************************ 

SUBROUTINE  DERIVat i ve ( KS I , N , I , XI , X2 , XDOT , MACH ) 
IMPLICIT  REAL*8 (A-H,0-Z) 
IMPLICIT  REAL*8(M) 
IMPLICIT  REAL*8 (K) 

DIMENSION  XDOT (N) , XI (N) , X2 (N) , KSI (N) , MACH (N) 
GO  TO  (20,20,10, 10, 10)  ,1 
10   XI (I) =  KSI (I) **2+ (l.ODO-MACH(I) **2) **2 
X2 (I) -XI (I) ** (1.0D0/3.0D0) 
XDOT (I) =1.0D0/X2 (I) 
GOTO  3  0 

2  0   XI (I) =- (KSI (I) **2- (l.ODO-MACH(I) **2) **2) 

X2 (I) =-Xl (I) ** (1.0D0/3.0D0) 
XDOT(I) =1.0D0/X2 (I) 

3  0   RETURN 

END 
c* *****************  ***************************************** 

C  SUBROUTINE  NEWTONCOTES  * 

c* *********************  ************************************* 

SUBROUTINE  NEWTONCOTES ( KS I , XDOT , X , SAVEX , SAVED , H , N , INDEX , I ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

IMPLICIT  REAL*8 (M) 

IMPLICIT  REAL*8 (K) 

DIMENSION  X (N) , XDOT (N) , SAVEX (N) , SAVED (N) , KSI (N) 

GOTO  (1,2,3,4,5)  , index 

1  SAVEX ( I ) =X ( I ) 
SAVED ( I ) =XDOT ( I ) 

KSI (I) -KSI (I) +0.2  5D0*H 
RETURN 

2  SAVED(I)=7.0D0*SAVED(I) +32 . 0D0*XDOT ( I ) 
KSI (I) =KSI (I) +0.2  5D0*H 

RETURN 

3  SAVED(I)=SAVED(I) +12 . 0D0*XDOT ( I ) 
KSI(I)=KSI (I) +0.2  5D0*H 

RETURN 

4  SAVED (I) =SAVED(I) +32 . 0D0*XDOT ( I) 
KSI (I) =KSI (I) +0.25D0*H 

RETURN 

5  X(I)=SAVEX(I) +(H/9  0.0D0) * (SAVED (I) +7 . 0D0*XDOT ( I) ) 
RETURN 

END 
c*****  **********  ******************************************** 
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************************************************************ 

THIS  FORTRAN  PROGRAM  USES  NUMERICAL  INTEGRATION  BASED  ON 
THE  NEWTON-COTES  FORMULA  TO  SOLVE  NUMERICALLY  EQNS .  41  AND  43. 

THE  OUTPUT  IS  IN  THE  FORM  OF  ij  vs.  Y 
************************************************************ 

program  ETTA 
PARAMETER  (N=5) 
PARAMETER  (K=2000) 
IMPLICIT  REAL*8 (A-H,0-Z) 
IMPLICIT  REAL*8 (M) 
INTEGER  INDEX, I, J 

DIMENSION  Y(N) ,YDOT(N) , SAVEY(N) , SAVED (N) ,TF(N) ,MACH(N) 
DIMENSION  ETA(N) , ETA11 (K) ,Y11 (K) , Yl (N) , Y2 (N) 
MACH(l) =0.8D0 
MACH(2) =0.9D0 
MACHO)  =1.0D0 
MACH(4) =1.1D0 
MACH(5) =1.2D0 
H=0.01D0 
q*********************************************************** 

C  INTEGRAL  FINAL  LIMITS  * 

q*********************************************************** 

DO  7  1=1,5 

IF(MACH(I) .EQ.10D0)  THEN 

TF(I) =0.0D0 

GOTO  7 
END  IF 
IF (MACH ( I ) . GT . 1 . 0D0 )  THEN 

TF(I) = (- (l.ODO-MACH(I) **2) ) ** (1 . 0D0/3 . 0D0) 

GOTO  7 
END  IF 

TF(I) = (l.ODO-MACH(I) **2) ** (1 . 0D0/3 . 0D0) 
7    END  DO 

Q*********************************************************** 

C  MAIN  PROGRAM  * 

q*  *********************************************************  * 

DO  100  1=1,5 
J=l 

ETA(I) =14.0D0 

Y(I) =2.0D0/ (ETA (I) ** ( 1 . 0D0/2 . 0D0 ) ) 
10   IF(ETA(I) -H.GT.TF(I) )  THEN 
2  0     do  50  index=l,5 

CALL  DERI Vat ive ( ETA , N , I , Yl , Y2 , YDOT , MACH ) 
CALL  NEWTONCOTES ( ETA , YDOT , Y , S AVEY , SAVED , H , N , INDEX , I ) 
50    end  do 

ETA11 (J)=ETA(I) 
Y11(J)=Y(I) 
J=J+1 
GOTO  10 
END  IF 
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H=0.001D0 

IF (ETA (I) -H.GT.TF(I) )  THEN 
GOTO  2  0 

END  IF 
80    DO  90  L=J-1, 1, -1 

WRITE(20+i, *)  Yll (L) , ETA11 (L) 
90    END  DO 
100  END  DO 

STOP 

END 
c**  ************  ********************************************* 

c  SUBROUTINE  DERIVative  * 

c*  ********  ************************************************** 

SUBROUTINE  DERIVa t  ive ( ETA , N , I , Yl , Y2 , YDOT , MACH ) 
IMPLICIT  REAL*8 (A-H,0-Z) 
IMPLICIT  REAL*8 (M) 

DIMENSION  YDOT(N) , Yl (N) , Y2 (N) ,ETA(N) ,MACH(N) 
GO  TO  (20,20, 10,10,10)  ,1 
10   Yl (I) =ETA(I) **3+ (l.ODO-MACH(I) **2) 
Y2 (I) -Yl (I) ** (1.0D0/2.0D0) 
YDOT(I) =1.0D0/Y2 (I) 
GOTO  3  0 

2  0   Yl (I)=ETA(I) **3- (l.ODO-MACH(I) **2) 

Y2 (I)=Y1 (I)  ** (1.0D0/2.0D0) 
YDOT(I)=1.0D0/Y2 (I) 

3  0   RETURN 

END 
c**  *********  ************************************************ 

C  SUBROUTINE  NEWTONCOTES  * 

Q*  *********  *   ************************************************ 

SUBROUTINE  NEWTONCOTES ( ETA , YDOT , Y , S AVEY , SAVED , H , N , INDEX , I ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

IMPLICIT  REAL*8(M) 

DIMENSION  Y (N) , YDOT (N) , SAVEY (N) , SAVED (N) , ETA (N) 

GO  TO  (1,2,3,4,5)  , index 

1  SAVEY ( I ) =Y ( I ) 
SAVED(I)=YDOT(I) 
ETA(I)=ETA(I) -0.2  5D0*H 
RETURN 

2  SAVED(I)=7.0D0*SAVED(I) +32 . 0D0*YDOT ( I) 
ETA(I)=ETA(I) -0.25D0*H 

RETURN 

3  SAVED(I)=SAVED(I) +12 . 0D0*YDOT ( I ) 
ETA (I) =ETA(I) -0.25D0*H 

RETURN 

4  SAVED(I)=SAVED(I) +32 . 0D0*YDOT ( I) 
ETA(I)=ETA(I) -0.2  5D0*H 

RETURN 

5  Y(I)=SAVEY(I) + (H/9  0.0D0) * (SAVED (I) +7 . 0D0*YDOT (I) ) 
RETURN 

END 
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APPENDIX   B 
FORM  OF  Ol   AND  0 

The  expressions  for  a  and  0  as  introduced  in  [REF.  2]  are 

a  =  ±C1(1-Ml)2  (117) 

P  =  +C2  (1-M„2)  (118) 

for  0.8SM.S1.2  . 

Upon  inserting  eqns .  117  and  118  in  eqns .  16  and  21  yields 

I 
X  =  f-3 — (119) 

o  jz2  +  (i-Mi)2 

y~"   J      /         i=  (120) 

a  2  i/f/3  +  (i-Ad) 

-d-wf)  3  v 

These  relations  are  numerically  integrated  and  plotted  in 
Figs.  16  and  17.  We  observe  the  same  solutions  for  MBsl.O  as 
expected,  but  distinguish  a  different  set  of  solutions  for  the 
subsonic  case  Moo<1.0  . 

Applying  the  same  procedure  as  in  Chap.  IV,  we  easily 
compute  the  expression  for  the  boundary  surfaces 

dY  /2\2    1     F  ,y 

?  dx  (121) 


l/  fj3  +  (1-Ml)  l3/  M„2(Y  +  D 
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The  numerical  integration  of  eqn.  121  is  shown  in  Fig.  18  for 
various  M,,  (e.g.,  0.8,0.9,0.97,1.0,1.1,1.2).  The  diverging 
contours  of  the  subsonic  solutions  appear  to  give 
contradictory  results  once  investigating  the  pressure 
coefficient  behavior 


C  -    ~2 


p     a£(y+i: 


fj  J  lz   +  (1-M„2)2  +  (1-M„2)J       (122) 


along  the  boundary  surface.  This  implies  as  shown  in  Fig.  19 
that  the  boundary  surface  within  the  portion  of  Cp>1.0  is  not 
a  realistic  solution. 
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APPENDIX   C 
FIGURES 
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Figure  1.  Numerical  integration  of  eqns.  40  and  42. 
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Figure  2.  Numerical  integration  of  eqns.  41  and  4  3 


E(Y)  vs.  Y  for  M  =0.8-1.2 
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Figure  3.  Numerical  solution  of  E(Y). 


Figure  4.  Equation  71. 
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